!>\brief
!!
!! Representation of the harmonic map state..
!!
!! \n
!!
!! The harmonic map is represented by the director and the Lagrange
!! multiplyer.
!<----------------------------------------------------------------------
module mod_hmap_state

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_octave_io, only: &
   mod_octave_io_initialized, &
   write_octave

 use mod_state_vars, only: &
   mod_state_vars_initialized, &
   c_stv

 use mod_grid, only: &
   mod_grid_initialized, &
   t_grid

 use mod_base, only: &
   mod_base_initialized, &
   t_base

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_hmap_state_constructor, &
   mod_hmap_state_destructor,  &
   mod_hmap_state_initialized, &
   t_hmap_state, clear, &
   new_hmap_state,      &
   write_octave,        &
   source

 private

!-----------------------------------------------------------------------

! Module types and parameters

 ! public members

 !> State information for the harmonic map
 !!
 !! Some pointers to the finite element space description are also
 !! introduced; such description is shared by all the state
 !! variables.
 type, extends(c_stv) :: t_hmap_state
  type(t_grid), pointer :: grid => null()
  type(t_base), pointer :: base => null()
  !> The director representation is organized on two indexes: the
  !! first one for the (vector) finite element basis, the second one
  !! for the element.
  real(wp), allocatable :: d(:,:)
  !> The Lagrange multiplier is a scalar function on each element.
  real(wp), allocatable :: q(:)
 contains
  procedure, pass(x) :: incr
  procedure, pass(x) :: tims
  !!> This subroutine is overidden both because of efficiency and
  !!! because many compilers leak memory when working with polymorphic
  !!! temporaries
  !procedure, pass(x) :: inlt => f_inlt
  !procedure, pass(z) :: alt  => f_alt
  !procedure, pass(z) :: mlt  => f_mlt
  procedure, pass(z) :: copy
  procedure, pass(x) :: source
  procedure, pass(x) :: source_vect
 end type t_hmap_state

 ! private members

 interface clear
   module procedure clear_hmap_state
 end interface

! Module variables

 ! public members
 logical, protected ::               &
   mod_hmap_state_initialized = .false.
 ! private members
 character(len=*), parameter :: &
   this_mod_name = 'mod_hmap_state'

 interface write_octave
   module procedure write_hmap_state
 end interface write_octave

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_hmap_state_constructor()
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_kinds_initialized.eqv..false.) .or. &
    (mod_messages_initialized.eqv..false.) .or. &
   (mod_octave_io_initialized.eqv..false.) .or. &
  (mod_state_vars_initialized.eqv..false.) .or. &
        (mod_grid_initialized.eqv..false.) .or. &
        (mod_base_initialized.eqv..false.) ) then  
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_hmap_state_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   mod_hmap_state_initialized = .true.
 end subroutine mod_hmap_state_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_hmap_state_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_hmap_state_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_hmap_state_initialized = .false.
 end subroutine mod_hmap_state_destructor

!-----------------------------------------------------------------------
 
 subroutine new_hmap_state(obj,grid,base)
  type(t_grid), intent(in), target :: grid
  type(t_base), intent(in), target :: base
  type(t_hmap_state), intent(out) :: obj

   obj%grid => grid
   obj%base => base
   !allocate( obj%d(base%mk,grid%ne) )
   !allocate( obj%q(        grid%ne) )
   allocate( obj%d(2,1) )
   allocate( obj%q(  1) )

 end subroutine new_hmap_state
 
!-----------------------------------------------------------------------
 
 pure subroutine clear_hmap_state(obj)
  type(t_hmap_state), intent(out) :: obj

   nullify( obj%grid )
   nullify( obj%base )
   ! allocatable components are implicitly deallocated
 
 end subroutine clear_hmap_state
 
!-----------------------------------------------------------------------

 subroutine incr(x,y)
  class(c_stv),        intent(in)    :: y
  class(t_hmap_state), intent(inout) :: x

   select type(y); type is(t_hmap_state)

   x%d = x%d + y%d
   x%q = x%q + y%q

   end select
 end subroutine incr

!-----------------------------------------------------------------------

 subroutine tims(x,r)
  real(wp),            intent(in)    :: r
  class(t_hmap_state), intent(inout) :: x

   x%d = r*x%d
   x%q = r*x%q

 end subroutine tims

!-----------------------------------------------------------------------

! subroutine f_inlt(x,r,y)
!  real(wp),     intent(in) :: r
!  class(c_stv), intent(in) :: y
!  class(t_f_state), intent(inout) :: x
!
!   select type(y); type is(t_f_state)
!
!   x%f = x%f + r*y%f
!
!   end select
! end subroutine f_inlt
!
!!-----------------------------------------------------------------------
!
! !> Overloaded for efficiency
! subroutine f_alt(z,x,r,y)
!  real(wp),     intent(in) :: r
!  class(c_stv), intent(in) :: x, y
!  class(t_f_state), intent(inout) :: z
!
!   select type(x); type is(t_f_state)
!    select type(y); type is(t_f_state)
!
!   z%f = x%f + r*y%f
!
!    end select
!   end select
! end subroutine f_alt
!
!!-----------------------------------------------------------------------
!
! !> Overloaded for efficiency
! subroutine f_mlt(z,r,x)
!  real(wp),     intent(in) :: r
!  class(c_stv), intent(in) :: x
!  class(t_f_state), intent(inout) :: z
!
!   select type(x); type is(t_f_state)
!
!   z%f = r*x%f
!
!   end select
! end subroutine f_mlt
!
!!-----------------------------------------------------------------------

 subroutine copy(z,x)
  class(c_stv),        intent(in)    :: x
  class(t_hmap_state), intent(inout) :: z

   select type(x); type is(t_hmap_state)

   z%d = x%d
   z%q = x%q

   end select
 end subroutine copy

!-----------------------------------------------------------------------

 subroutine source(y,x)
  class(t_hmap_state), intent(in) :: x
  class(c_stv), allocatable, intent(out) :: y

   allocate(t_hmap_state::y)
 
   select type(y); type is(t_hmap_state)

   call new_hmap_state( y , x%grid , x%base )

   end select
 end subroutine source
 subroutine source_vect(y,x,m)
  integer, intent(in) :: m
  class(t_hmap_state), intent(in)  :: x
  class(c_stv), allocatable, intent(out) :: y(:)

  integer :: i

   allocate(t_hmap_state::y(m))

   select type(y); type is(t_hmap_state)

   do i=1,m
     call new_hmap_state( y(i) , x%grid , x%base )
   enddo

   end select
 end subroutine source_vect

!-----------------------------------------------------------------------

 subroutine write_hmap_state(f,var_name,fu)
  integer, intent(in) :: fu
  type(t_hmap_state), intent(in) :: f
  character(len=*), intent(in) :: var_name
 
  character(len=*), parameter :: &
    this_sub_name = 'write_f_state'

   write(fu,'(a,a)')    '# name: ',var_name
   write(fu,'(a)')      '# type: struct'
   write(fu,'(a)')      '# length: 2' ! number of fields

   ! field 01 : d
   write(fu,'(a)')      '# name: d'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(f%d,'<cell-element>',fu)

   ! field 02 : q
   write(fu,'(a)')      '# name: q'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(f%q,'c','<cell-element>',fu)

 end subroutine write_hmap_state

!-----------------------------------------------------------------------

end module mod_hmap_state

